import numpy as np
import matplotlib.pyplot as plt
def run():
    # 任意磁导率下的单层最优吸波介电常数计算
    fmax = 18  # 频率最大值  单位 Ghz
    fmin = 2  # 频率最小值 单位 Ghz
    d = 3  # 给定的任意边长厚度 单位mm
    e1max = 500  # 相对介电常数实部最大值默认1500，如遇到2GHz下算的值相对偏大，可按照趋势修改
    e1min = 0
    e2max = 100  # 相对介电常数虚部，默认100，基本不需要修改
    e2min = 0
    u1 = 1  # 相对磁导率实部为1
    u2 = 0  # 相对磁导率虚部为0
    c = 299792458  # 真实真空光速，采用这个
    c1 = 3e8  # 真空近似光速
    df = 0.08  # 频率步长，按照实际要求修改，可以更小
    de1 = 1  # 介电实部步长，无需修改
    de2 = 1  # 介电虚部步长，无需修改
    # 初始值定义
    f = np.arange(fmin, fmax + df, df)  # 频率节点 单位Ghz
    e1 = np.arange(e1min, e1max + de1, de1)  # 介电实部取值
    e2 = np.arange(e2min, e2max + de2, de2)  # 介电虚部取值
    z = np.zeros((len(e1), len(e2)), dtype=complex)
    RL = np.zeros((len(e1), len(e2)))
    RLmin = np.zeros((len(f), 1))
    row = np.zeros_like(RLmin)
    col = np.zeros_like(RLmin)
    RL1 = np.zeros_like(RLmin)
    e1m = np.zeros_like(RLmin)
    e2m = np.zeros_like(RLmin)
    e11 = np.zeros((201, 1))
    e21 = np.zeros_like(e11)
    z1 = np.zeros((201, 201), dtype=complex)
    smallRL1 = np.zeros((201, 201))
    e12 = np.zeros_like(e11)
    e22 = np.zeros_like(e11)
    z2 = np.zeros_like(z1)
    smallRL2 = np.zeros((201, 201))
    row1 = np.zeros((len(f), 1))
    col1 = np.zeros((len(f), 1))
    row2 = np.zeros((len(f), 1))
    col2 = np.zeros((len(f), 1))
    [er1, er2] = np.meshgrid(e1, e2)
    for k in range(0, len(f)):
        z = np.sqrt((u1 - u2 * 1j) / (er1 - er2 * 1j)) * np.tanh(
            1j * 2 * d * 1e-3 * np.pi * f[k] * 1e9 / c * np.sqrt((er1 - er2 * 1j) * (u1 - u2 * 1j)))
        RL = 20 * np.log10(np.abs((z - 1) / (z + 1))).T
        if k == 0:
            [row[k], col[k]] = np.where(RL == np.nanmin(RL))
        elif e1[int(row[k - 1])] + 1 > e1max:
            [row[k], col[k]] = np.where(RL == np.nanmin(RL))
        else:
            RL1 = RL[0:int(row[k - 1]) + 2, :]
            [row[k], col[k]] = np.where(RL1 == np.nanmin(RL1))
        # 精确度到相对介电常数小数点后2位
        e11 = np.arange(e1[int(row[k])] - de1, e1[int(row[k])] + de1 + de1 / 100, de1 / 100)
        e21 = np.arange(e2[int(col[k])] - de2, e2[int(col[k])] + de2 + de2 / 100, de2 / 100)
        [er11, er21] = np.meshgrid(e11, e21)
        z1 = np.sqrt((u1 - u2 * 1j) / (er11 - er21 * 1j)) * np.tanh(
            1j * 2 * d * 1e-3 * np.pi * f[k] * 1e9 / c * np.sqrt((er11 - er21 * 1j) * (u1 - u2 * 1j)))
        smallRL1 = 20 * np.log10(np.abs((z1 - 1) / (z1 + 1))).T
        [row1[k], col1[k]] = np.where(smallRL1 == np.nanmin(smallRL1))
        # 精确度到相对介电常数小数点后4位
        e12 = np.arange(e11[int(row1[k])] - de1 / 100, e11[int(row1[k])] + de1 / 100 + de1 / 10000, de1 / 10000)
        e22 = np.arange(e21[int(col1[k])] - de2 / 100, e21[int(col1[k])] + de2 / 100 + de2 / 10000, de2 / 10000)
        [er12, er22] = np.meshgrid(e12, e22)
        z2 = np.sqrt((u1 - u2 * 1j) / (er12 - er22 * 1j)) * np.tanh(
            1j * 2 * d * 1e-3 * np.pi * f[k] * 1e9 / c * np.sqrt((er12 - er22 * 1j) * (u1 - u2 * 1j)))
        smallRL2 = 20 * np.log10(np.abs((z2 - 1) / (z2 + 1))).T
        [row2[k], col2[k]] = np.where(smallRL2 == np.nanmin(smallRL2))
        e1m[k] = e12[int(row2[k])]
        e2m[k] = e22[int(col2[k])]
        RLmin[k]=np.nanmin(smallRL2)

    plt.figure()
    plt.plot(f, e1m)
    plt.plot(f, e2m)
    plt.xlabel('Frequency /GHz')
    plt.ylabel('er')
    plt.show()

    plt.figure()
    plt.plot(f, RLmin)
    plt.xlabel('Frequency /GHz')
    plt.ylabel('Reflection Loss/ dB')
    plt.show()

if __name__=='__main__':
    run()